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Abstract 

This paper presents a general formulation of the CIP/multi-moment finite volume method (CIP/MM FVM) for arbi- 
trary order of accuracy. Reconstruction up to arbitrary order can be built on single cell by adding extra derivative 
moments at the cell boundary. The volume integrated average (VIA) is updated via a flux-form finite volume formula- 
tion, whereas the point-based derivative moments are computed as local derivative Riemann problems by either direct 
interpolation or approximate Riemann solvers. 
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1. Introduction 

The multi-moment concept underlying the CIP method (Cubic-Interpolated Pseudo-particle or Constrained Inter- 
polation Pro file) 1 16] provides a general methodology to construct numerical schemes with great flexibility. One of 
the major outcome from the practice so far to implement the multi-moments in computational fluid dynamics is that 
we can build high order schemes on a relatively compact grid stencil using multi-moments, and these moments can 
be carried forward in time separately by completely different numerical approaches. 

Some schemes have been developed for practical use based on VIA and SIA (Surface-Integrated Average) lfT3ll lif 
and on VIA and PV (Point Value) yDlj]. The later is much more suitable for unstructured or other complex com- 
putational grids where a point-wise local Riemann problem can be posed at any specified point to update the PV. It 
is found that increasing the number of the PVs is a simple way to get higher order schemes. We have devised and 
verified the schemes up to 4th order on 2D triangular unstructured grid for both scalar and system conservation laws 
by employing both VIA and PV moments. On the other hand, making use of the first derivative at the cell boundary 
as another moment has been ever used in the so-called CIP-CSL4(CIP-Conservative Semi-Lagrangian with 4th order 
polynomial) advection scheme @]. 

We in this paper explore further the possibility to construct conservative CIP/multi-moment formulation of arbi- 
trary order over single cell using more derivative moments. The spatial reconstruction based on multi-moments is 
described in section 2. The numerical formulation for scalar hyperbolic conservation law is presented in section 3. 
The extension to Euler equations is discussed in section 4. Section 5 ends the paper with a few conclusion remarks. 

2. The multi-moment spatial reconstruction 

The essential point in high resolution scheme is how to reconstruct the interpolation function to find the numerical 



flux at the boundary of each grid cell. Among the most widely used are, for example, the MUSCL scheme [12], 
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Dynamics, ISBN 978-4-339-04597-0, publish by Coronasha, 2009. 
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the ENO scheme[2] and the WENO scheme [5]. In all of these schemes the interpolation is based only on the cell- 
averaged values of the physical field to be reconstructed. In this section, we describe a numerical interpolation that 
makes use of not only the volume-integrated average over each mesh cell but also the derivatives at the cell boundary. 
We call the present formulation the "multi-moment" reconstruction to distinguish it from the aforementioned ones 
which should be more properly refer to as the "single-moment" reconstruction. 

We consider a physical field variable <p{x) over a one-dimensional domain divided into control volumes (mesh 
cells) [x,-i/ 2 ,*, + i/2]; i = 1,2,...,/. 

There is a flexibility in choosing the discretised moments for field variable <p(x, t). The primary moment of the 
finite volume method is the volume-integrated average (VIA) over each mesh cell 



Ax, J x , 



(x, t)dx (1) 



where Ax,- = x, + i/2 - x,-i/2- The spatial derivatives up to Kt\\ order at cell boundary 

^ /+ . = §(%j,f)^ft (2) 
with fc = 0, 1, - - • , K; 

are also used as the moments. Note that the point value (PV) P 4>i+in is actually equivalent to the 0th derivative moment 
D *(°)0 i+ 1 . We will denote the kth derivative of <p in respect to any variable d k <f>/d/3 k , by di ^<p occasionally hereafter. 

Given one VIA v (p i and 2(K + 1) derivative moments Djt ®0 i± i over [*;-i/2, x /+ i/2], as well as the first-order deriva- 
tive or gradient di that is computed in terms of other independent moments, we can construct a [2(K + 1) + l]th order 
cell-wise polynomial €>,(x, t) with 2(K + 2) constrained conditions as follows, 



X' 



+I <D,(x, t)dx = v (p h (3) 



df 4>i(Jt i+ i) = D *™(p i+ i, k = 0,1, ••• ,K; (4) 



df 3>i(*,--i) = D <«0,._! , k = 0, 1, •• • ,K; (5) 
d^i(Xi) = dt. (6) 
Thus, the piecewise interpolation polynomial, 

K 

Q> i {x) = Y J ak{x-x i _i % ) k , k = Q,l,--,K (7) 

jfc=0 

is constructed over cell i. All the coefficients can be uniquely computed from (O,©,© and ((6]). The first-order 
derivative or gradient of the interpolation function di can be approximated in terms of the known moments. For 
example, a [2(K + 1) + l]th order polynomial is obtained if we specify di - d x Q>i(Xi) with 4>,(jc) is computed from 
constraint conditions ([3), (|4]i and (0. Furthermore, a slope limiting can be imposed to d, to suppress the numerical 
oscillation (see [ 15] for details). We used a single-cell minmod limiter in this paper. 

It is obvious that the reconstruction discussed above can have a 2(K + l)th order accuracy for smooth solutions. 

3. The scalar conservation laws 

In this section, we consider the scalar conservative law as follows, 

dcf> , 8f{cp) 



dt dx 



= 0, (8) 



where <p is the scalar state variable and f(<p) is the flux function. Assuming the hyperbolicity, we have a real charac- 
teristic velocity, u = d^f(<p). 





Figure 1: Advection transport of a harmonic wave after 400 steps (CFL=().l) by schemes with derivative moments of 0th order (top-left), 1th order 
(top-right), 2nd order (bottom-left) and 3rd order (bottom-right). 



The governing equations for the derivative moments can be directly derived from (O as, 



a«(^ i+ .) = -af«(/(^4 +r 

k = 0, ,K; 



(9) 



where / is the numerical flux function consistent to /. 

It is observed from the reconstruction that the derivatives moments up to the Kth order D *<£>0 and flux 3, / 
are continuous. Thus, we can update the derivative moments D '^(p i+ i for k = 0, 1, • • • ,K — 1 by (TTOb with the 

spatial derivatives of the flux function d% f i+ 1 directly computed from the derivative moments that readily defined and 
computed at the cell interface as 



(10) 



When one advances the highest order derivative moment D > W<^> , d <K+1) f / dx K+l is required, which, however, might 
not be continuous at the cell boundaries. We make use the the simple Lax-Friedrichs splitting in terms of the spatial 
derivatives of the flux function and the state variable as 



(id 



where a i+ 1 is the largest value of the characteristic speed in the related region. 
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Figure 2: Advection transport of a square pulse after 2000 steps (CFL=0.1) by schemes with derivative moments of Oth order (top-left), 1th order 
(top-right), 2nd order (bottom-left) and 3rd order (bottom-right). 



The state variables (p^ , and <f)P , are computed from the multi-moment reconstructions separately built for 
cells [jc,_i /2 ,x, + i /2 ] and [x !+1/2 , x, +3/2 ], i.e. 

4>f +k = <S>i(x i+ i) and = 
The corresponding (X + 1 ) derivatives are 

<9f +1) 0f + i = af +1) ^i + p and (12) 
e + Xi=d? +1) * !+ i(*, + i). (13) 



It should be noted that we have used an assumption similar to ||9]| jUlli ] in getting a homogeneous and linearized 
Riemann problem for spatial derivatives for the state variable. 

In order to update the VIA moment, we integrate © over [*;-i/2, *,-+i/ 2 ], yielding the following conservative 
formulation, 

= _±. in mj i+ £ - f ( ^ i+ ,)), (H) 

where f denotes the numerical flux at cell boundary, and is computed directly from the derivative moments readily 
updated at the cell boundaries. 

The semi-discretized time evolution equations (9) and (TBI are predicted in time by a TVD|01 or a 4th order 
Runge-Kutta method. At every substep, we first update the derivative moments by (9), and then use these updated 
moments to evaluate the flux function in (TT4l . 

We computed an advected harmonic wave with a wavelength of 20/7Ax. Fig.l shows the results for reconstruc- 
tions using derivatives moments up to different orders. It is found that even a short wave can be adequately resolve if 
higher order derivative moments are used. 



Fig. 2 reveals the effect of limiting. The numerical oscillation associating discontinuities are eliminated by the 
slope switching that is also constructed within a single mesh element. 



4. The Euler conservation laws 

In this section, the numerical formulation presented above is implemented to the inviscid Euler conservation laws. 
The conservative form of the one-dimensional Euler equations is written as follows, 



at 



<w 

dx 





f 9 1 
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0, u = 


pu 


, F = 


pu 2 + p 




, e , 




, u(e + p) , 



(15) 



where U is the vector of conservative variables and F is the vector of inviscid fluxes. Denoted by p is the density, 
u the velocity, e the total energy and p the pressure that is obtained by the equation of state for the perfect gas 
p = (e - pu 2 1 '2){y - 1). The ratio of the specific heats y is specified as 1.4 in this paper. 
The volume-integrated average (VIA) moment over mesh cell i, 



AX; J x 



U(x, t)dx, 



(16) 



and the derivative moments up to Kt\\ order at cell boundary, 



DAk) U i+ l = -gj( X H>V With * = 0,1,- 



,K; 



(17) 



are treated as the model variables. 
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Figure 3: Numerical results of Sod's problem at t = 0.2 separately computed by the multi-moment reconstructions with K = 0, 1, 2 and the 5th-order 
WENO scheme are plotted. 



Analogous to the scalar case, the governing equations for the derivative moments are derived from (TT~5T > as, 



5%pWU i+ i)= (18) 
_ d (k+D ( F (o*(*),U)). + , , k = 0, 1, • • ■ , K; 
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where F i s the ve ctor of the numerical flux function consistent to F. 

Given D '*'U i+ i continuous for k = 0, 1, ■ • • ,K,we can update the derivative moments D '®U i+ i for k = 0, 1, ■ • • , K— 

1 by ( fT9l ) with the spatial derivatives of the flux function d x F evaluated directly from 



5®F( fl *C0)u i+i ,o«(i)U i+ i,- • • , D >»U i+ i). 



(19) 



Similar to the scalar case, 8% might be not continuous at the cell boundaries. So, the highest order derivative 
moment D *®U has to be solved from a Riemann problem in terms of the spatial derivative. To this end, we use the 
linearization assumption in j^ lfioll and write the (K + l)th spatial derivative of the flux function as 



df F = Adf U, 



(20) 



where A is the Jacobian matrix. 

The conventional flux splitting algorithms can be adopted here in terms of the spatial derivative quantities. 

The (K + l)th derivatives of the state variables and the flux function at the cell interface are computed from the 
multi-moment reconstructions separately built over two neighboring cells as (12) and (13) component-wisely in terms 
of the the state variables or characteristic variables. 
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Figure 4: The close-up views for shock (left) and contact discontinuity(right). 



The VIAs of the conservative variables U on cell i are updated by integrating ( TT~5b over [JCf-i/i, Xf+1/2], which 
results in a finite volume formulation, 

Given the derivative moments D »®U i+ ! at cell boundaries, the numerical fluxes ¥( Dx ^V i+ i) in the above equation are 
directly found. 

Again, the Runge-Kutta method is used for time integration for all moments. 

A ID shock tube test[7] was computed to verify the present method for Euler conservation laws. We include the 
numerical result of the 5th-order WENO scheme[5] as well for comparison. Shown in Fig. 3, the numerical results of 
the present scheme with different orders are quite competitive. The close-up plots for shock and contact discontinuity 
are given in Fig. 4. Both linear and non-linear discontinuities are well resolved with correct locations. It is observed 
that better resolution can be obtained by simply increasing the order of the derivative moments. 
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5. Concluding remarks 



A formulation that uses high order derivative moments has been suggested and tested. Given all the derivative 
moments that are continuous at cell boundaries and updated separately, the resulting numerical formulation is still 
single-cell based and quite computationally efficient. In case that the derivative moments are defined and continuous 
at the cell boundary, the numerical fluxes can be computed directly as in the IDO scheme [1], while for the spatial 
derivative higher than the continuous one, we simplify and cast it into a linearized derivative Riemann problem^. 
The present formulation is substantially different from the ADER methodic lioll 111] where all the derivatives are 
discontinuous at cell boundaries, thus is more efficient. 

Our numerical results show that the resolution of the scheme can be improved by simply increasing the order 
of the derivative moments involved. With the simple slope limiting |l5|, the numerical oscillation around the large 
gradient can be effectively suppressed. 

Although the multi-dimension implementation remains an open problem to be further explored, one can expect 
the present scheme as an accurate efficient solver for ID conservation laws. 
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